n-ta odmocnina a isqrt
Otázka od: Tomas Krysl
28. 11. 2002 21:10
Jeste jednou zdravim,
tohle je tak trochu OT. Potrebuji realizovat n-tou odmocninu celeho cisla
(vysledek ma byt tez cele cislo) bez pouziti FP aritmetiky, tj. bez
logaritmu a exp (procesor HC12 - tj. rychlost je vsechno). Nekde jsem mel
napsany algoritmus a nemuzu ho najit. Jednalo se o jakysi genialni
rekurentni vzorecek. Pokud se Vam neco takoveho vali v kompu, prosim poslete
(na prog. jazyku nezalezi).
Dal jsem v nejakych C-ckovskych snipetech nasel nasledujici vypocet
2.odmocniny (sqrt) a moc mi to nerika. Pokud jste se tim nekdo zabyval,
prosim poslete mi vysvetleni (pry jsem se to mel ucit na zakladce, ale asi
jsem zrovna chybel). Nerad bych ve svem programu mel neco, cemu nerozumim.
Dik.
T.
/* +++Date last modified: 05-Jul-1997 */
#include <string.h>
struct int_sqrt {
unsigned sqrt,
frac;
};
#define BITSPERLONG 32
#define TOP2BITS(x) ((x & (3L << (BITSPERLONG-2))) >> (BITSPERLONG-2))
/* usqrt:
ENTRY x: unsigned long
EXIT returns floor(sqrt(x) * pow(2, BITSPERLONG/2))
Since the square root never uses more than half the bits
of the input, we use the other half of the bits to contain
extra bits of precision after the binary point.
EXAMPLE
suppose BITSPERLONG = 32
then usqrt(144) = 786432 = 12 * 65536
usqrt(32) = 370727 = 5.66 * 65536
NOTES
(1) change BITSPERLONG to BITSPERLONG/2 if you do not want
the answer scaled. Indeed, if you want n bits of
precision after the binary point, use BITSPERLONG/2+n.
The code assumes that BITSPERLONG is even.
(2) This is really better off being written in assembly.
The line marked below is really a "arithmetic shift left"
on the double-long value with r in the upper half
and x in the lower half. This operation is typically
expressible in only one or two assembly instructions.
(3) Unrolling this loop is probably not a bad idea.
ALGORITHM
The calculations are the base-two analogue of the square
root algorithm we all learned in grammar school. Since
we're
in base 2, there is only one nontrivial trial multiplier.
Notice that absolutely no multiplications or divisions are
performed.
This means it'll
be fast on a wide range of processors.
*/
void usqrt(unsigned long x, struct int_sqrt *q)
{
unsigned long a = 0L; /* accumulator */
unsigned long r = 0L; /* remainder */
unsigned long e = 0L; /* trial product */
int i;
for (i = 0; i < BITSPERLONG; i++) /* NOTE 1 */
{
r = (r << 2) + TOP2BITS(x); x <<= 2; /* NOTE 2 */
a <<= 1;
e = (a << 1) + 1;
if (r >= e)
{
r -= e;
a++;
}
}
memcpy(q, &a, sizeof(long));
}
#ifdef TEST
#include <stdio.h>
#include <stdlib.h>
main(void)
{
int i;
unsigned long l = 0x3fed0169L;
struct int_sqrt q;
for (i = 0; i < 101; ++i)
{
usqrt(i, &q);
printf("sqrt(%3d) = %2d, remainder = %2d\n",
i, q.sqrt, q.frac);
}
usqrt(l, &q);
printf("\nsqrt(%lX) = %X, remainder = %X\n", l, q.sqrt, q.frac);
return 0;
}
#endif /* TEST */
-----
Ing. Tomáš Krýsl,
programátor a analytik
ALVA Strakonice s.r.o.
Písecká 893, 386 01 Strakonice
++420 377 430 701
GSM: 776 892 206
tomkrysl@quick.cz
Odpovedá: Petr Brant
29. 11. 2002 8:55
Da se pouzit napriklad Newtonova metoda tecen (ale pozor, selhava v
blizkosti cisla 0, pak je vhodnejsi pouzit jinou metodu, viz napr.
http://pascal.fjfi.cvut.cz/%7Elimpouch/numet/). Podrobnejsi popis algoritmu
(format doc) posilam na soukromy mail.
RNDr. Petr Brant [brant@dcomm.cz]
http://web.redbox.cz/petr.brant
D&COMM s.r.o.
Korunovační 6
Praha 7
tel. +420724007234
Subject: n-ta odmocnina a isqrt
tohle je tak trochu OT. Potrebuji realizovat n-tou odmocninu celeho cisla
(vysledek ma byt tez cele cislo) bez pouziti FP aritmetiky, tj. bez
logaritmu a exp (procesor HC12 - tj. rychlost je vsechno). Nekde jsem mel
napsany algoritmus a nemuzu ho najit. Jednalo se o jakysi genialni
rekurentni vzorecek. Pokud se Vam neco takoveho vali v kompu, prosim poslete
(na prog. jazyku nezalezi).